# Erik Fast, Quan Lam, Zhuoyu Li, Kyle Tiffany
# Math 381 Project, Fall 2011

simulate <- function(bag.employees, nonbag.employees)
{
  # the total number of fans that enter through each gate
  # order: south, southwest, west, northwest, north, student
  gate.totals = c(3763, 2811, 12587, 9593, 5864, 2034)
  
  # total fans to enter the stadium, according to the data
  # this includes other gates that we are not modeling
  total.fans <- 53644
  
  # get the proportion of people that pass through each gate
  gate.proportions <- gate.totals / total.fans
  
  # the number of fans employees can admit during one time interval
  entries.per.bag.employee <- 21.7
  entries.per.nonbag.employee <- 63.1
  
  gate.flow.bag <- bag.employees * entries.per.bag.employee 
  gate.flow.nonbag<- nonbag.employees * entries.per.nonbag.employee
  
  # expected number of total fans for each of the time intervals
  #Case 1: Average attendance
  #total.expected <- c(1227, 1487, 1550, 2007, 2351, 2602, 3136, 3634, 3950, 4290, 4398, 4137, 3676, 2839, 2030)
  #Case 2: Full Attendance
  total.expected <- c(1894, 2222, 2077, 2772, 3267, 3771, 4456, 5014, 5129, 5175, 4600, 3887, 2664, 1718, 1186)
  #Case 3: Low attendance
  #total.expected <- c(812, 1069, 1122, 1509, 1746, 1943, 2426, 3046, 3320, 4057, 4389, 4481, 4496, 4498, 3875, 2914)
  # keep track of the number of people waiting at a gate
  residual <- rep(0, 6)
  
  for (expected in total.expected)
  {
    # generate a random number of arrivals from Poisson distribution
    total.arrivals <- rpois(n=1, lambda=expected)
    
    # distribute the arrivals among the gates based on the proportions 
    gate.arrivals <- gate.proportions * total.arrivals + residual
    
    # calculate the number of fans that are not let in during the interval
    a1<-gate.arrivals*0.9 - gate.flow.bag
    a2<-gate.arrivals*0.1 - gate.flow.nonbag
    for (i in 1:6) {
    	if (a1[i]<0) a1[i]<-0
    	if (a2[i]<0) a2[i]<-0
    }

    residual <- a1+a2

    # we don't want negative residuals
    residual <- c(max(0, residual[1]), max(0, residual[2]), max(0, residual[3]), max(0, residual[4]), max(0, residual[5]), max(0, residual[6]))
  }
  
  return(residual)
}

#Original staff assignment
#base.bag.employees <- c(10, 8, 20, 12, 19, 10)
#base.nonbag.employees <- c(7, 6, 22, 10, 10, 8)

#Improved staff assignment
base.bag.employees <- c(15, 12, 37, 19, 21, 16)
base.nonbag.employees <- c(2, 2, 5, 3, 3, 2)

#adjust the number of staff at all gates
#case 1: Average attendance
#improved.bag.employees <- base.bag.employees + c(-6,-6,-8,+3,-8,-8)
#improved.nonbag.employees <- base.nonbag.employees + c(-1,-1,-1,0,-1,-1)
#Total -38

#case 2: Full attendance
improved.bag.employees <- base.bag.employees + c(-4,-3,-6,+4,-6,-5)
improved.nonbag.employees <- base.nonbag.employees + c(-1,-1,-1,0,-1,-1)
#Total -25

#case 3: Low attendance
#improved.bag.employees <- base.bag.employees + c(-4,-5,-4,+6,-6,-6)
#improved.nonbag.employees <- base.nonbag.employees + c(-1,-1,-1,0,-1,-1)
#Total -24

num.simulations <- 10000
base.residual <- rep(0, 6)
improved.residual <- rep(0, 6)

for (i in 1:num.simulations)
{
   base.residual <- base.residual + simulate(base.bag.employees, base.nonbag.employees)
   improved.residual <- improved.residual + simulate(improved.bag.employees, improved.nonbag.employees)
}

base.residual <- base.residual / num.simulations
improved.residual <- improved.residual / num.simulations

base.residual
improved.residual